static char help[] = "This program calculates the 1D hard-sphere (HS) chemical potential as \
defined by the Density Functional Theory (DFT) developed by Rosenfeld.\n\n";
#include <dft.h>

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc, char *argv[])
{
  MPI_Comm       comm;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = PetscInitialize(&argc, &argv, (char *) 0, help);CHKERRQ(ierr);
  comm = PETSC_COMM_WORLD;

/*   PetscInt  dim; */
/*   PetscReal L[3]; */
/*   PetscReal h[3]; */
  
  DA da;
/*   ierr = DACreate3d(comm, DA_NONPERIODIC, DA_STENCIL_STAR,  */
/*                     10,10,10, */
/*                     PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE, */
/*                     NUM_SPECIES, 0, */
/*                     PETSC_NULL, PETSC_NULL, PETSC_NULL, */
/*                     &da);  */
/*   CHKERRQ(ierr); */
  ierr = DACreate1d(comm, DA_NONPERIODIC, 
                    10,
                    NUM_SPECIES, 0,
                    PETSC_NULL, 
                    &da); 
  CHKERRQ(ierr);
  ierr = DASetUniformCoordinates(da, -5.0, 5.0, 
                                 PETSC_NULL, PETSC_NULL, 
                                 PETSC_NULL, PETSC_NULL); 
  CHKERRQ(ierr);
  PetscReal R[NUM_SPECIES];
  for(unsigned int i = 0; i < NUM_SPECIES; ++i) {
    R[i] = 1.0;
  }
  PetscReal T = 273.0;
  PetscReal z[NUM_SPECIES];
  DFT *dft;
  ierr = DFTCreate(da, R, z, T, &dft); CHKERRQ(ierr);

  Vec rho, muHS;
  ierr = DAGetGlobalVector(da, &rho); CHKERRQ(ierr);
  ierr = DAGetGlobalVector(da, &muHS); CHKERRQ(ierr);
  ierr = VecSet(rho, 0.1); CHKERRQ(ierr);
  ierr = DFTCalculateMuHS(dft, rho, muHS); CHKERRQ(ierr);
  ierr = DARestoreGlobalVector(da, &rho); CHKERRQ(ierr);
  ierr = DARestoreGlobalVector(da, &muHS); CHKERRQ(ierr);

  ierr = PetscFinalize(); CHKERRQ(ierr);
  PetscFunctionReturn(0);
}// main()
